Self-replication and splitting of domain patterns in reaction-diffusion systems with 

fast inhibitor 
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An asymptotic equation of motion for the pattern interface in the domain-forming reaction- 
diffusion systems is derived. The free boundary problem is reduced to the universal equation of 
non-local contour dynamics in two dimensions in the parameter region where a pattern is not far from 
the points of the transverse instabilities of its walls. The contour dynamics is studied numerically 
for the reaction-diffusion system of the FitzHugh-Nagumo type. It is shown that in the asymptotic 
limit the transverse instability of the localized domains leads to their splitting and formation of the 
multidomain pattern rather than fingering and formation of the labyrinthine pattern. 
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I. INTRODUCTION 

Pattern formation is a remarkable phenomenon typi- 
cal for many physical, chemical, and biological systems 
inside and outside of thermal equilibrium As a 

rule, these are extremely complicated phenomena. How- 
ever, in many cases pattern formation may be explained 
on the basis of systems of reaction-diffusion equations 
of the activator-inhibitor type. The systems of this 
kind include electron-hole and gas plasma, semiconductor 
and superconductor structures, systems with uniformly 
generated combustion material, chemical reactions with 
auto-catalysis and cross-catalysis, models of morphogen- 
esis and population dynamics (see, for example, [^]-[7j and 
references therein). The simplest example of such a sys- 
tem is a pair of reaction-diffusion equations: 



f)f) 
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where 9 is the activator, r\ is the inhibitor, I and L are the 
characteristic length scales, and tq and are the char- 
acteristic time scales of the activator and the inhibitor, 
respectively, q and Q are certain non-linear functions, 
and A is the bifurcation parameter. 

Kerner and Osipov showed that the properties of pat- 
terns and pattern formation scenarios in the systems de- 
scribed by Eqs. (0) and (0) are chiefly determined by 
the parameters e = l/L and a = rg/r^ and the shape 
of the nullcline of the equation for the activator 
In many cases this nullcline is N-shaped (Fig. 1). In 
such N systems static domain patterns may form when 
Kl HHf§- These patterns are essentially the domains 
of high and low values of the activator separated by the 
narrow walls whose width is of order I. Recent experi- 
ments and numerical simulations have revealed a lot of 
new pattern formation scenarios in these systems such 
as growth of fingers, tip splitting, spot replication, and 



formation of labyrinthine patterns [jL0-15|. These effects 
are associated with the fact that at certain parameters 
the patterns become unstable with respect to the trans- 
verse perturbations. Muratov and Osipov developed a 
general asymptotic theory of such instabilities in an ar- 
bitrary N system described by Eqs. (§) and @ §]. They 
have shown that the instabilities are determined by the 
motion of the pattern walls. Goldstein, Muraki, and Pet- 
rich derived an equation of motion for a simple reaction- 
diffusion system of FitzHugh-Nagumo type in the limit of 
fast inhibitor and small activator-inhibitor coupling and 
drew similar conclusions about the transverse instabili- 
ties of the domain patterns |H. The numerical simula- 
tions of concrete reaction-diffusion systems in the limit 
of fast inhibitor (a ^> e) showed that the instabilities 
typically lead to the formation of labyrinthine patterns. 
Spot replication was observed only in the case of slow 
inhibitor (a < e) 

Numerical solution of Eqs. (||) and (§) for a concrete 
model shows that for the same values of the parameters 
there exist qualitatively different types of stable static 
solutions. These are the solutions in the form of the 
labyrinthine patterns [Fig. 2(a)] and the multidomain 
patterns [Fig. 2(b)]. The labyrinthine pattern forms as 
a result of the instability of a single domain which is 
taken as an initial condition in Fig. 2(a). In the sim- 
ulation of Fig. 2(b) the initial condition was taken in 
the form of a random arrangement of domains. As a 
result, a metastable multidomain pattern forms in the 
system. We emphasize that the only difference between 
the simulations of Fig. 2 is the initial conditions, all pa- 
rameters characterizing the system itself are identical in 
both cases. 

The aim of this paper is to show that in the asymptotic 
limit e — > the transverse instabilities of the pattern 
walls will always lead to splitting of domains and for- 
mation of the multidomain patterns in the limit of fast 
inhibitor. To show this we will reduce the partial differ- 
ential equations problem in an arbitrary d-dimensional N 
system to the free boundary problem in the limit e — > 
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FIG. 1. Nullclines of Eqs. (1) and (2) with q = 6 3 - 6 - r) and Q = 9 + 77 - A. 



for arbitrary a. We will further reduce this problem near 
the instability points of the domain patterns to the prob- 
lem of non-local contour dynamics for a e and obtain 
a universal equation of motion for the interface in two 
dimensions. This equation will be studied numerically 
for a concrete model. 



II. FREE BOUNDARY PROBLEM 

In this section it is convenient to measure the lengths 
in the units of I and time in the units of Te, respectively. 
Then Eqs. (g) and (g) become: 

— = A9-q(0,ri,A), (3) 

a- 1 ^ = e- 2 Ar,-Q(e,r,,A). (4) 

Ohta, Mimura, and Kobayashi developed an approach 
which allowed them to reduce equations similar to Eqs. 



(g) and (g) to the problem of the interface dynamics in 
the case of slow inhibitor (a < e) in the limit e — > |fl6| . 
Here we will use their approach to derive the equations 
of the interface dynamics for Eqs. (g) and (g) and show 
that these equations are essentially the same in the case 
of fast and slow inhibitor. 

Let us introduce the local coordinate system in the 
vicinity of the interfaces of the pattern (Fig. g). For 
a point x let p be the distance from the point x to the 
interface, and the (d — l)-dimensional coordinate £ the 
projection of x on the submanifold p — const, p is as- 
sumed to be positive if the point is in the region of high 
activator values and negative otherwise. The distribu- 
tion of the activator varies on the length scale 1 near 
the interface. Since the characteristic length of the vari- 
ation of the inhibitor is (T 1 ^> 1, 77 can be considered 
constant in the direction perpendicular to the interface: 
V = F r0lrL the general considerations follows that 

the curvature of the interface K(^) has to be small [g,g. 



Therefore, to the first power in e and K we may write 
Eq. (g) as follows: 

The sign of K is such that it is positive if the interface is 
convex into the low activator region. Notice that £ enters 
only as a parameter in this equation. For times t ^> 1 
any solution of Eq. (g) connecting high and low values 
of the activator will become close to the solution in the 
form of a front propagating with the velocity v in the 
direction perpendicular to the interface. For any given 
point £ one can then write 

,. s d9 d 2 9 r r,^d9 .„ , ,v, 



where we introduced the self-similar variable z = —p — vt. 
The velocity is positive when it points out of the high ac- 
tivator value domain. The boundary conditions for this 
equation are: 

0(+oo)=0ii(O, 0(-oo)=0 i3 (O: (7) 

where 9n^ are the minimal and the maximal roots of the 
equation 

q(0*AOMO) = (8) 

for any given £. The third solution 0*2 (£) 01 Eq. (g) 
which lies between On and 9a may be used to fix the po- 
sition of the interface relative to 9(z) by requiring that 
9 = ^i2(0 on the interface. Thus, represents the 
normal velocity of the pattern interface at a point £. 
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FIG. 2. Two types of extended domain patterns: labyrinthine pattern (a) and multidomain pattern (b). Results of the numer- 
ical simulations of Eqs. (1) and (2) with q — O 3, — 8 — 1} and Q — 8 + rj — A. The parameters used are: e = 0.05, a — 0.2, A — —0.4. 
The system size is 20L x 20L. The time of the simulation is 500t,,. For explanations of the initial conditions see the text. 



The velocity v can be found by solving Eq. (g) with the 
abovementioned boundary conditions. The consistency 
condition which is obtained from Eq. (||) by multiply- 
ing it by d9/dz and integrating over z gives the following 
expression for the velocity in terms of 9(z): 



v = -K- 



JJlq[%m)de 

J—oo KdzJ 



(9) 



From this equation follows that the front velocity v is of 
order 1. Notice that Eq. (||) is essentially an equation of 
motion of a particle in the potential 



Uo 



q(0,Vi)dO 



(10) 



with the friction proportional to v + K with time z |}] ^| . 
In N systems Ug is a double hump potential (for those 
values of 77, for which it has only a single hump Eq. (^) 
has no solutions), so Eq. (||) has a unique solution which 
satisfies the boundary conditions in Eq. (0) for the given 
values of r\i and K, which corresponds to the trajectory 
going from the top of one hump to the top of the other. 
Thus, the velocity of the interface is a single- valued func- 
tion of the curvature and the value of 77 at the interface. 

Far from the interfaces the activator varies on the char- 
acteristic length e -1 , so the Laplacian in Eq. (||) in that 
region can be dropped and Eq. (j^) reads: 



I") 



So, the relation between 9 and 77 far from the interfaces 
is local in space. Since 77 varies on the length scales of 
order e _1 or if -1 and the velocity of the front is of order 



1, the characteristic time of the variation of the interface 
velocity is also e -1 > 1 or K -1 3> 1, what justifies the 
use of Eq. (j|) for finding the distribution of the activator 
near the interfaces. The characteristic time of variation 
of the inhibitor field caused by the interface motion will 
also be e _1 or K^ 1 in the case a ^> e (fast inhibitor), 
or a^ 1 , which is even greater than the previous in the 
case a < e (slow inhibitor), so the time derivative in Eq. 
( |TT| ) is less or of order e and, therefore, can be neglected. 
Then, the values of 9 and 77 far from the interfaces are 
simply related by the equation of local coupling: 



#,r?)=0. 



(12) 



In view of Eq. (|lj), the inhibitor must satisfy the equa- 
tion: 



a- 1 ^ = e- 2 A V -Q(9( V ),r 1 ), 



(13) 



with the following boundary conditions at the interfaces 
(written in terms of the local curvilinear coordinates p 
and £): 



(14) 



The dependence 6(rj) in Eq. (|l3| ) is the solution of Eq. 
(|l2l). This dependence is multi- valued, so for the re- 
gions of high and low activator values one should take the 
branches connected with 9^ and 9n, respectively. Thus, 
we recovered the result obtained by Ohta, Mimura, and 
Kobayashi in Ref. j[6| in the case of fast inhibitor as well. 
The term in the left-hand side of Eq. ( p"3| ) is of order e/a 
which can be set equal to zero in the case of fast inhibitor. 

Equation ( |l3| ) is essentially nonlinear since it involves 
the inverse function of essentially nonlinear function 
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FIG. 3. The schematics of the system's geometry (two dimensions). The thick solid line shows the interface. e p and e$ are 
the local orthogonal basis of the curvilinear coordinate system p, £; p indicates the distance from a given point to the interface, 
the dashed line indicates the surface p = const. 



q(9,rj), even if the Eq. (0) were linear. Moreover, the 
right-hand side of Eq. (pjf) becomes singular at the 
points #o and 9' (see Fig. [j]) where q'g(9o,r]o) — and 
q'g(9' ,r]' ) = 0, respectively, since the function rj(9) be- 
comes singular at these points. When the distribution 
of 1] at a point xq reaches the values of t)q or ri' , a sud- 
den down-jump or up-jump (local breakdown), occurs at 
that point, respectively 15|. In terms of the free 

boundary problem this corresponds to the instantaneous 
creation of a new interface at the point Xq, which will 
start to evolve according to the equation of the inter- 
face motion. This is an important ingredient of the free 
boundary problem describing the dynamics of the pat- 
tern interface that should not be left out in solving this 
problem. Notice that local breakdown is responsible for 
the domain splitting in the one-dimensional N systems 
11- 

Equations (§) - (|) together with Eqs. (Q - (0) and 
the rule for dealing with the singularities of Eq. ( |lli| ) 
discussed in the previous paragraph are the closed set of 
equations which defines completely the dynamics of any 
domain pattern in the system described by Eqs. (^) and 
(|J) in the limit e — > both in the case of slow and fast 
inhibitor. These equations were obtained with the ac- 
curacy to e < 1 and K <C 1. Note that the presence 
of the curvature term in the equation of motion for the 
interface suggests only that the interface has to be suffi- 
ciently smooth, more exactly, it has to be smooth on the 
length scales of order 1. This, of course, does not mean 
that the interfaces are not allowed to intersect, fuse, or 
loose their connectivity in the course of their evolution. 
Self-replicating domains, which will be studied in the sub- 
sequent sections, are a good example of the latter. Also, 
a domain may disappear if it shrinks too much. 



III. NON-LOCAL CONTOUR DYNAMICS 

As was shown in Ref. ||, when a ^> e static domain 
patterns may be stable only when their characteristic size 
is much smaller than the characteristic length of the in- 
hibitor variation because of the stabilizing action of the 
interaction between the walls of the pattern. Typically, 
a pattern destabilizes when the characteristic distance 
between its walls (in the units of the previous Section) 
becomes of order e~ 2 / 3 ||. In this situation the veloc- 
ity of the non-stationary pattern interface must be much 
smaller than one. Also, the value of the inhibitor must 
be close to the value of rj s at which the velocity of the 
pattern wall is equal to zero in one dimension. According 
to Eq. (||) with K — 0, the value of r] s must satisfy 



q(9,r, s ,A)d9 = Q, q(9 sh3 ,r, s ,A) = 0, (15) 



where 9 S ± and 9 S ^ are some constants, which generally 
depend on A. This situation takes place when the bi- 
furcation parameter A is close to the values of Ab or A' b 
where the characteristic distance between the walls of the 
domains with high or low values of the activator becomes 
zero in the limit e — » 0. 

For definiteness let us consider a single domain of high 
values of the activator. Then the value of the activator 
inside the domain will be close to 9 s s, and to 9 S \ outside. 
Let us introduce the new variables: 

9 = 9-9 sl , fj = v-Vs (16) 

Equations (|l2|) and (|l|) can then be linearized. In the 
case of fast inhibitor we will have: 



1 e (9 s i,r) s )9 + q'J9 s i,r) s )fj = 0, 



(17) 



= e- 2 Afi-Cfj-aI{x)-Q{9 sl , Vs ), (18) 



where 



C = Q' v (6 sl ,r ]s ) 



Qe(8si,ris)q' n (0si,r]s) 



(19) 



(20) 



and I(x) is the indicator function which is equal to 1 
if x is inside the domain and outside. The value of 
Q(9 s i l ru) is small for A close to A&. Note that in writing 
Eq. ( |l8| ) we neglected the piecewise-constant potential 
inside the domains which is present upon the lineariza- 
tion in the general case since one has to linearize Eq. ( |l2] ) 
inside and outside of the domain separately. However, 
this is justified when the domain size is much smaller 
than the characteristic length of the inhibitor variation 
since the potential then is only a perturbation and can 
be neglected in the zeroth order. 

Since the velocity of the front is small, it can also be 
linearized in f). This, however, extremely simplifies the 
problem since one no longer has to solve the non-linear 
eigenvalue problem in Eq. (|^). Indeed, in the linear ap- 
proximation in f) we can replace 9{z) in the denominator 
of Eq. (||) by the solution in the form of the static one- 
dimensional front and expand the numerator in 77, so we 
immediately obtain for the velocity at a point £ on the 
interface: 

v(0 = -K(0 - UOZ- 1 [" q' n {6, Vs)dd, (21) 



and 



Z= / ^2(U e (9 sl ,7 ls )-Ue(0,Vs))de, (22) 



where Ue is defined in Eq. (|10|). 

The distribution of fj(x) at each moment in time and, 
therefore, the values of f] on the interface, which deter- 
mine the interface velocity, can be found by solving Eq. 
( |l8| ) by means of the Green's function: 

fj( x ) = _ a J G{ x _ x ')l( x ')d d x '. (23) 

Specifically, in the infinite two-dimensional system 



G{x - x') = ^-K Qx - x'\eVC), 
2n 



(24) 



where Kq is the modified Bcsscl function. Following the 
idea of Ref. jlJJ , let us transform the integral in Eq. j2^ ) 
in two dimensions into the contour integral along the 
interface. Using the defining equation for the Green's 
function and the Gauss theorem, we find: 



fj(x) 



where r is the vector from the point x to x' , where the 
point x' lies on the interface (see Fig. |||), and the inte- 
gration is over the interface, K\ is the modified Bessel 
function. In writing Eq. ( |25| ) it was taken into account 
that the surface integral obtained from Eq. ( pi| ) can as 
well be written in terms of the vector product involving 
the tangent vector dr. The interface has to be oriented 
counterclockwise in order to get the sign right. If there is 
more than one domain in the system, one has to add up 
the contributions of each domain given by the integral in 

For a single domain of the size of order e -2 / 3 one can 
expand the Bessel function in Eq. (pq), so we obtain the 
following equation of motion for the interface: 



v(0 = -K(0- 



where 



BQ{d sl ,Vs) 



aCZ 



Be 2 
~A^Z 



(ln(0.54e VC) + In r)r x dr, (26) 



B = -(Q(8 s3 ,r is )-Q(6 sl ,i ls )) q' v (e, Vs )d9, (27) 

Jo sl 

and the point x is now on the interface. If we introduce 
the renormalized coordinates, and time and introduce the 
new control parameter A: 



x = xe 2 '\ t = te i '\ A=- 



BQ(6 sl ,r] s )e 
aCZ 



-2/3 



(28) 



we will eliminate the e-dependence in Eq. (|2^) (except 
for the weak logarithmic dependence). Thus, we obtained 
the asymptotic equation of motion for a localized do- 
main with the characteristic size much smaller than the 
characteristic size of the inhibitor variation. This equa- 
tion was derived with the accuracy to e 2 / 3 . Note that at 
this point all information about the specific nonlincaritics 
of the system is contained only in a few numerical con- 
stants: B, C and Z. These constants are all positive for 
the reaction-diffusion systems of the activator-inhibitor 
type Furthermore, they can be incorporated into the 
rescaled e = eZC 3 ^ 2 /B, provided that time and coordi- 
nate are also suitably rescaled. So, the dynamics of the 
localized domains in any two-dimensional N system not 
far from Ai, depends only on A and e, the last dependence 
being logarythmically weak. 

Equation (^) is also good for the description of several 
interacting domains if the distance between them is not 
much greater than e -2 / 3 . In order to study the interac- 
tion of domains separated by the distances of order e _1 
one has to use Eqs. (21) and (p5|). 

Equation ( p6[ ) contains a large logarithm which mul- 
tiplies the integral § r x dr which up to a coefficient is 
the area of the domain. This means that with the log- 
arithmic accuracy the area of the evolving domain has 
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to be conserved. Also, the non-local term in Eq. ( |26| ) is 
an increasing function of r, so the parts of the interface 
which are close to each other attract, and the parts which 
are far from each other repel. Then the instability which 
appears for a radially-symmetric domain when its radius 
reaches certain value |9| cannot lead to the growth of a 
labyrinthine pattern, but rather will result in the domain 
splitting. The domains will then go apart and the pro- 
cess of splitting will repeat, until the systems becomes 
filled with the multidomain pattern. This mechanism of 
the instability development is qualitatively different from 
the one discussed in Ref. IH. 



IV. SELF-REPLICATING DOMAINS IN A 
CONCRETE REACTION-DIFFUSION SYSTEM 

In this Section we will apply the results obtained above 
and show that the proposed mechanism of the instability 
development indeed takes place in a concrete reaction- 
diffusion system. To do this we will use the system which 
is described by Eqs. (|J) and (||) with 



Q = 9 + r]-A. 



(29) 



(30) 



The nullclines of this system for A = —0.6 are shown in 
Fig. 1. 

This system is particularly convenient because the de- 
pendences of its characteristic parameters on A are es- 
pecially simple. Also, the non-linear eigenvalue problem 
in Eq. (^) in this system admits exact solution. Simple 
calculation shows that: 

T) s = 0, 0.1 = -1, 0,3-1, (31) 
and the coefficients involved in Eq. ( p6[ ) are: 

o = 2, B = 4, C=| Z=^. (32) 



The value of Q(6 s i,r) s ) becomes zero at A = A\, with 
Ai, = —1. Also, the homogeneous state of the system 



0,, 



-\A\ 



1/3 



Vh 



-|^4| 1 / 3 (1 - \A\ 2 / 3 ) is stable for 



A < A c = -1/3^3. 

The solution of Eq. (||) can be sought in the form 
6{z) = atanh&z + c, where a, 6, and c are constants. As 
a result, the dependence of the front velocity on rji (with 
K = 0) is implicitly given by the following equation: 



v ir 

9 



(33) 



< •v/3/2, the maximum value 
T)o , where rio = 2/3 V3- Thus, 



The velocity v satisfies 
being achieved at r\ = r/o, wnere rip 
for all values of r\ for which Eq. (y_3|) is not singular the 
dependence v{r\i) is single- valued. 

Another special feature of this system is the fact that 
in it the piecewise-constant potential mentioned in Sec. 



Ill is identically zero, so Eqs. (El]) and (25) also describe 



the dynamics of any complex pattern with the character- 
istic domain size much smaller than (T 1 for any values 
of A. These include the complex patterns forming in the 
late stages of the development of the transverse instabili- 
ties P|l5|]. In fact, the analog of Eq. fl25|), non-local both 
in space and time, may be obtained also in the case of ar- 
bitrary a. This equation, together with Eq. (pl|), could 
be used to study the pulsations (breathing) of complex 
domain patterns in this system. 

The direct numerical simulations of Eqs. @ and (^) 
is a formidable task. The main difficulty here is the fact 
that for small e there are two very different length scales, 
so the simulations require enormous amounts of compu- 
tational power. Recently, it became possible to perform 
extensive numerical simulations of the system under con- 
sideration using massive parallelization on a supercom- 
puter jl5| . The authors were able to simulate the system 
of sufficiently large sizes with e ~ 0.05. The values of 
e ~ 0.01 are already very hard to simulate even on the 
very fast computer. 



The main result of these simulations relevant to the 
present discussion is that as a result of the development of 
the transverse instability a localized domain transforms 
into a labyrinthine pattern, which is all connected if A is 
not far from Af,. This effect takes place at e > 0.01. For 
those values of e domain splitting and self-replication was 
observed only when a < e | fL5f . We performed numeri- 
cal simulations of this system for e < 0.01 and saw that 
a destabilizing localized domain actually splits into two. 
However, the simulation becomes excessively long at this 
point, so it is impossible to see if the forming domains 
will split in turn. 

Equation (p6|), from the other hand, allows to study 
the interfacial dynamics for arbitrarily small e. We per- 
formed numerical simulations of Eq. (Eq) for the values of 



the parameters when a single localized domain becomes 
unstable with respect to transverse perturbations of its 
walls. 

Figure ^ shows the evolution of the almost circular 
domain when its radius is close to the critical radius 
at which the domain loses stability with respect to the 
m = 2 mode. Initially the domain elongates, but at 
some point the distance between the walls becomes so 
small that the domain splits into the three disconnected 
pieces. The resulting domains continue to grow until the 
larger domains split again into seven (not shown in the 
Figure), and the process goes on. Notice that in order 
for this process to take place, the value of e has to be 
very small. As is seen from the numerical simulations of 
Eqs. (||) and (|j), when e is not very small the domain 
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FIG. 4. Destabilization of the circular domain. Results of the numerical simulation of Eq. (26) for the considered model 
with A = 20 and e = 10 -4 . The box indicates the size of 30 x 30. The length and time are measured in the rescaled units given 
by Eq. (28). 



can elongate and transform into a stripe without split- 
ting. This effect also takes place in the reaction-diffusion 
systems with the weak activator- inhibitor coupling pj . 

Figure |5| shows the results of another simulation when 
the system is further away from the instability point. 
There the initial stage of the domain growth is similar to 
the one shown in Fig. ^. However, the walls in the center 
come closer to each other than in the previous case, so 
at some point the domain splits into two. The resulting 
domains continue to grow and split in turn. This process 
is essentially self-replication of the domains, as a result of 
which the system will become filled with the multidomain 
pattern. 

We emphasize that unless the system is very close to 
the point A = Af,, the localized domains will always be 
unstable and, once excited, will transform into a multido- 
main pattern via self-replication even in the case of fast 
inhibitor, if e is small enough. This is a completely uni- 
versal result independent of any property of the system. 
In particular, this conclusion does not depend on whether 
the system is monostable or bistable. As was noted in the 
previous Section, the only nontrivial system-dependence 
is contained in the logarithmic term, and, therefore, this 
dependence is weak. 

The results obtained for the system under considera- 
tion do not change in a wide range of e. When e > 0.01, 



Eq. (26) seizes to be a good approximation for Eqs. pi] ) 
and (25). For extremely small e (e < 10~ 9 ) there exist a 
narrow region of the values of A close to the critical value 
at which the disk transforms into a stable domain in the 
form of a dumbbell. When the value of A is increased, 
the domain self-replication will occur even for such small 
values of e. 



V. CONCLUSION 

Thus, we have shown that in the case of fast inhibitor 
the transverse instability of the localized domains in the 
reaction-diffusion systems with N-shaped nullcline for the 
activator will always lead to domain splitting and forma- 
tion of multidomain pattern rather than the formation of 
labyrinthine patterns, if e is sufficiently small. 

This effect was in fact observed in a quasi two- 
dimensional experiment with the current filaments form- 
ing in n-GaAs in the process of avalanche breakdown JL7) . 
There the radially-symmetric current filaments destabi- 
lized and split as the current in the sample increased and 
the radius of a filament grew. At some critical value 
of current a filament split into two and the filaments 
that formed split in turn until their radii become suffi- 
ciently small. Non-symmetrically distorted or elongated 
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FIG. 5. Self-replication of the domains. Results of the numerical simulation of Eq. (26) for the considered model with 
A = 30 and e = 10~ 4 . The box indicates the size of 30 x 30. The length and time are measured in the rescaled units given by 
Eq. (28). 



filaments were not observed. 

Hagberg and Meron suggested that splitting of do- 
mains is the consequence of the parity-breaking front 
transitions [noncquilibrium Ising-Bloch (NIB) transi- 
tions] associated with the variations of the curvature of 
the domain walls Jl3| . However, their arguments may 
apply only to the bistable reaction-diffusion systems in 
which the inhibitor is slow. We showed here that the 
splitting of domains in fact occurs in the systems with 
the fast inhibitor and is determined by the non-local in- 
teraction between different portions of the domain inter- 
face, regardless of whether the system is monostable or 
bistable, so as a rule the NIB transitions should not be 
responsible for domain splitting. It is clear that in the 
systems with the slow inhibitor domain splitting will oc- 
cur even easier since the inhibitor will not be able to react 
on the motion of the walls of the domains locally not only 
in space, but also in time. This is also confirmed by the 
direct numerical simulations of Eqs. (0) and (^) fig] . 

Notice that the condition e -C 1 itself is the necessary 
condition for the existence of the static domain patterns 
in the considered systems |^-|e|, so, in fact, the transition 
from a localized domain to the multidomain pattern con- 
sisting of disconnected localized domains filling up the 
space must be the major mechanism of the transverse 
instability development. Multidomain patterns were in- 
deed observed in the chemical systems [n0| and in the 



high-frequency gas discharge Eg ]. Nevertheless, numeri- 
cal simulations and experimental observations also show 
that for small but finite e localized domains may trans- 
form into extended labyrinthine patterns |To|- |T^ , |l4| , |T5| . 
This does not seem to be the consequence of the finite 
width of the interface, but rather a peculiarity of the in- 
hibitor dynamics. This effect can actually be controlled 
by varying the strength of interaction between the ac- 
tivator and the inhibitor [the constant B in Eq. (p6|)]. 
Goldstein, Muraki, and Petrich analyzed the equation 
similar to Eqs. ([H]) and ( |25| ) in the case B ~ e in the 
limit e — * and found that the transverse instability 
leads to the formation of the connected labyrinthine pat- 
tern. This suggests that by changing e or the coupling 
strength B one could control whether the multidomain 
pattern or the labyrinthine pattern will form as a result 
of the domain instability. 

It is important to note that the free boundary prob- 
lem obtained from Eqs. (|l|) and (||) in the limit e -c 1 
contains considerably less information about the nonlin- 
earities of the system than the initial partial differential 
equations problem. Moreover, according to Eq. (p6|), the 
behavior of any localized pattern not far from the points 
Af, or A' b is universal in the sense that the dynamics of 
the interface can be described by only a few renormal- 
ized parameters. This universality was discussed earlier 
in the context of the instabilities of the domain patterns 
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P). This suggests that the free boundary formulation 
of the pattern dynamics might be a more advantageous 
starting point for dealing with the problems of domain 
pattern formation rather than the formulation in terms 
of reaction-diffusion equations. 

The phenomenology of pattern formation similar to 
the one discussed in the present article is also observed 
in various equilibrium systems with competing interac- 
tions (see, for example, || and references therein). No- 
tice that the reaction-diffusion system described by Eqs. 
(|l|) and (Q) in the limit r n — ► describes the kinetics of 
a system with competing interactions, if q = f{&) — r) 
and Q = rj + B8, where B is a constant and f(9) is 
some cubic-like function [ p^| , pO|Jll| ] . These equations de- 
scribe, for example, the kinetics of microphase separa- 
tion of block copolymers . The universality of the 
results obtained above suggests that self-replication of 
domains must be a common feature of the systems with 
competing interactions in the case of repulsive long-range 
interactions of Coulombic type and strong separation of 
length scales. 

The author is grateful to V. V. Osipov for valuable 
discussions. 
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